library(matrixStats)
library(sp)
library(maptools)
## Checking rgeos availability: TRUE
library(maps)
library(raster)
library(colorRamps)
library(rgdal)
## rgdal: version: 1.2-6, (SVN revision 651)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-4
library(scales)
library(animation)
library(dismo)
library(deldir)
## deldir 0.1-14
library(rgeos)
## rgeos version: 0.3-23, (SVN revision 546)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.2-4
## Polygon checking: TRUE
library(gstat)
Here I transform the breeding timing data into a SpatialPointsDataFrame and then create a presence/absence matrix to show which birds are breeding on a given date.
load("Global_bird_breeding_dates_julian.Rdata")
# change dataframe into SpatialPointsDataFrame
spring_data <- data.frame(spring_data)
class(spring_data[,5]) <- "numeric"
class(spring_data[,6]) <- "numeric"
nonas <- spring_data[which(is.na(spring_data[,5])==FALSE | is.na(spring_data[,6])==FALSE),]
nonas <- nonas[which(is.na(nonas[,5])==FALSE & is.na(nonas[,6])==FALSE),]
coord <- coordinates((nonas[,6:5]))
bird_pts <- SpatialPointsDataFrame(data =nonas, coords = coord, proj4string = CRS("+proj=longlat +datum=WGS84"))
# make a presence absence matrix
namevector <- vector(length = 365)
for (i in 1:365){
namevector[i] <- paste0("day", i)
}
dates_matrix <- data.frame(matrix(ncol = 365, nrow = 1780))
dim(dates_matrix)
colnames(dates_matrix) <- namevector
dates_matrix[] <- 0
bird_pts_binded <- spCbind(bird_pts, dates_matrix)
bird_pts <- bird_pts_binded
i <- namevector[1]
for (i in namevector){
bird_pts[which(bird_pts$Egg_start_julian_date <= as.numeric(substr(i, 4, nchar(i))) & as.numeric(substr(i, 4, nchar(i))) <= bird_pts$Fledge_end_julian_date), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 1
bird_pts[which(bird_pts$Egg_start_julian_date > bird_pts$Fledge_end_julian_date), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 3
bird_pts[which(bird_pts[[i]] == 3 & as.numeric(substr(i, 4, nchar(i))) >= bird_pts$Egg_start_julian_date), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 1
bird_pts[which(bird_pts[[i]] == 3 & as.numeric(substr(i, 4, nchar(i))) <= bird_pts$Egg_end_julian_date), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 1
bird_pts[which(bird_pts[[i]] == 3), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 0
}
bird_pts$day1 <- 0
plot(bird_pts[which(bird_pts$day1 == 1),])
save("bird_pts", file = "bird_points_presence_absence.Rdata")
load("bird_points_presence_absence.Rdata")
Here I plot points representing breeding birds throughout the year on top of the 1941-1950 raster and on top of the 2007-2016 raster.
list_dates <- seq(as.Date("2017-01-01"), as.Date("2017-12-31"), by = "day")
list_dates <- format(list_dates, format = "%b %d")
list_files_south <- list.files("Past Interpolation Southern", full.names = TRUE)
list_files_north <- list.files("Past Interpolation Northern", full.names = TRUE)
for (i in 1:365){
list_files_south[i] <- paste0("Past Interpolation Southern/" , i, "past_northern.Rdata")
list_files_north[i] <- paste0("Past Interpolation Northern/" , i, "past_northern.Rdata")
}
list_files_south_mod <- list.files("Present Interpolation Southern", full.names = TRUE)
list_files_north_mod <- list.files("Present Interpolation Northern", full.names = TRUE)
for (i in 1:365){
list_files_south_mod[i] <- paste0("Present Interpolation Southern/" , i, "present_southern.Rdata")
list_files_north_mod[i] <- paste0("Present Interpolation Northern/" , i, "present_northern.Rdata")
}
c1 <- colorRampPalette(c("navy", "blue", "magenta2", "indianred1", "yellow", "white"), bias = 8.5)
saveGIF({
for (i in 1:365){
par(mfrow=c(1,2))
load(as.character(list_files_south[i]))
plot(masked, main = paste0("1941-1950 Growing Degree Day Accumulation | ", list_dates[i] ), mar = c(2, 4, 1, 2) +0.1, col = c1(5000), breaks = seq(0,5000, by = 100), legend = FALSE, xlim = c(-180, 180), ylim = c(-90,90), asp = 1)
abline(h = 0, lty = 2)
load(as.character(list_files_north[i]))
plot(masked, add = TRUE, col = c1(5000), breaks = seq(0,5000, by = 100), legend = FALSE)
plot(bird_pts[which(bird_pts[[paste0("day", i)]]==1),], add = TRUE, pch = 16, cex = 0.5, col = "lawngreen")
load(as.character(list_files_south_mod[i]))
plot(masked, main = paste0("2007-2016 Growing Degree Day Accumulation | ", list_dates[i] ), mar = c(2, 4, 1, 2) +0.1, col = c1(5000), breaks = seq(0,5000, by = 100), legend = FALSE, xlim = c(-180, 180), ylim = c(-90,90), asp = 1)
abline(h = 0, lty = 2)
load(as.character(list_files_north_mod[i]))
plot(masked, add = TRUE, col = c1(5000), breaks = seq(0,5000, by = 100), legend = FALSE)
plot(bird_pts[which(bird_pts[[paste0("day", i)]]==1),], add = TRUE, pch = 16, cex = 0.5, col = "lawngreen")
print(i)
}
}, movie.name = "global_interpolation_birds.gif", ani.width = 1000, ani.height = 500, interval = 0.3)